library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.4
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 3.6.2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Get bird fragmentation data

bird_frag <- read.csv(file = "https://mlammens.github.io/ENS-623-Research-Stats/data/Logan_Examples/Chapter9/Data/loyn.csv")
summary(bird_frag)
##      ABUND            AREA            YR.ISOL          DIST       
##  Min.   : 1.50   Min.   :   0.10   Min.   :1890   Min.   :  26.0  
##  1st Qu.:12.40   1st Qu.:   2.00   1st Qu.:1928   1st Qu.:  93.0  
##  Median :21.05   Median :   7.50   Median :1962   Median : 234.0  
##  Mean   :19.51   Mean   :  69.27   Mean   :1950   Mean   : 240.4  
##  3rd Qu.:28.30   3rd Qu.:  29.75   3rd Qu.:1966   3rd Qu.: 333.2  
##  Max.   :39.60   Max.   :1771.00   Max.   :1976   Max.   :1427.0  
##      LDIST            GRAZE            ALT       
##  Min.   :  26.0   Min.   :1.000   Min.   : 60.0  
##  1st Qu.: 158.2   1st Qu.:2.000   1st Qu.:120.0  
##  Median : 338.5   Median :3.000   Median :140.0  
##  Mean   : 733.3   Mean   :2.982   Mean   :146.2  
##  3rd Qu.: 913.8   3rd Qu.:4.000   3rd Qu.:182.5  
##  Max.   :4426.0   Max.   :5.000   Max.   :260.0

Log transform AREA

bird_frag$LOG_AREA <- log(bird_frag$AREA)
x_values <- bird_frag$LOG_AREA %>% 
  round(2)
y_values <- bird_frag$GRAZE %>% 
  round(2)
z_values <- bird_frag$ABUND %>% 
  round(2)
# Define regression plane -------------------------------------------------
# Construct x and y grid elements
x_grid <- seq(from = min(bird_frag$LOG_AREA), to = max(bird_frag$LOG_AREA), length = 50)
y_grid <- seq(from = min(bird_frag$GRAZE), to = max(bird_frag$GRAZE), length = 50)
# Construct z grid by computing
# 1) fitted beta coefficients
# 2) fitted values of outer product of x_grid and y_grid
# 3) extracting z_grid (matrix needs to be of specific dimensions)
beta_hat <- bird_frag %>% 
  #lm(ABUND ~ LOG_AREA + GRAZE, data = .) %>% 
  lm(ABUND ~ LOG_AREA + GRAZE + log10(DIST) + log10(LDIST) + YR.ISOL + ALT, data = .) %>%
  coef()

fitted_values <- crossing(y_grid, x_grid) %>% 
  #mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid)
  mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid + 
           beta_hat[4]*mean(log10(bird_frag$DIST)) + beta_hat[5]*mean(log10(bird_frag$LDIST)) +
           beta_hat[6]*mean(bird_frag$YR.ISOL) + beta_hat[7]*mean(bird_frag$ALT))

z_grid <- fitted_values %>% 
  pull(z_grid) %>%
  matrix(nrow = length(x_grid)) %>%
  t()
# Define text element for each point in plane
text_grid <- fitted_values %>% 
  pull(z_grid) %>%
  round(3) %>% 
  as.character() %>% 
  paste("ABUND: ", ., sep = "") %>% 
  matrix(nrow = length(x_grid)) %>%
  t()
# Plot using plotly -------------------------------------------------------
plot_ly() %>%
  # 3D scatterplot:
  add_markers(
    x = x_values,
    y = y_values,
    z = z_values,
    marker = list(size = 5),
    hoverinfo = 'text',
    text = ~paste(
      "ABUND: ", z_values, "<br>",
      "GRAZE: ", y_values, "<br>",
      "LOG_AREA: ", x_values 
    )
  ) %>%
  # Regression plane:
  add_surface(
    x = x_grid,
    y = y_grid,
    z = z_grid,
    hoverinfo = 'text',
    text = text_grid
  ) %>%
  # Axes labels and title:
  layout(
    title = "3D scatterplot and regression plane",
    scene = list(
      zaxis = list(title = "y: ABUND"),
      yaxis = list(title = "x2: GRAZE"),
      xaxis = list(title = "x1: LOG_AREA")
    )
  )